function H = generate_channel_matrix(Mx, My, d, P, K, angles)
% GENERATE_CHANNEL_MATRIX 生成3D信道矩阵（文献Section II）
%   输入：
%     Mx, My: URA尺寸（x/y方向天线数）
%     d: 天线间距（波长倍数）
%     P: 每条路径数
%     K: 总用户数
%     angles: 用户角度结构体（generate_user_angles输出）
%
%   输出：
%     H: [M×K]信道矩阵，M=Mx*My
%
%   文献依据：公式(1)-(3)

    % ========== 1. 参数验证 ==========
    validate_parameters(Mx, My, d, P, K, angles);
    
    % ========== 2. 初始化信道矩阵 ==========
    M = Mx * My; % 总天线数
    H = zeros(M, K); % 初始化信道矩阵
    
    % ========== 3. 遍历所有用户 ==========
    for k = 1:K
        % 3.1 获取当前用户角度信息
        gamma_x = angles(k).gamma_x; % [P×1]
        gamma_y = angles(k).gamma_y; % [P×1]
        
        % 3.2 生成路径增益（公式1）
        % 复高斯分布：实部/虚部独立N(0,1/√(2P))
        g = (randn(P,1) + 1i*randn(P,1)) * sqrt(1/(2*P));
        
        % 3.3 初始化用户信道向量
        h_k = zeros(M, 1);
        
        % ========== 4. 遍历所有路径 ==========
        for p = 1:P
            % 4.1 计算x方向导向矢量（公式3）
            phase_x = -2*pi*d*gamma_x(p)*(0:Mx-1)';
            e_x = exp(1i*phase_x) / sqrt(Mx);
            
            % 4.2 计算y方向导向矢量（公式3）
            phase_y = -2*pi*d*gamma_y(p)*(0:My-1)';
            e_y = exp(1i*phase_y) / sqrt(My);
            
            % 4.3 计算完整导向矢量（Kronecker积）
            e_xy = kron(e_x, e_y); % [M×1]
            
            % 4.4 累加路径贡献
            h_k = h_k + g(p) * e_xy;
        end
        
        % ========== 5. 存储用户信道 ==========
        H(:, k) = h_k;
    end
    
    % ========== 6. 信道归一化 ==========
    % 确保E[|h_k^{m,n}|^2]=1（文献Section II）
    H = H / sqrt(mean(abs(H(:)).^2));
end

%% ========== 子函数：参数验证 ==========
function validate_parameters(Mx, My, d, P, K, angles)
    % 天线参数
    if ~isscalar(Mx) || Mx < 1 || mod(Mx,1) ~= 0
        error('Mx must be a positive integer');
    end
    if ~isscalar(My) || My < 1 || mod(My,1) ~= 0
        error('My must be a positive integer');
    end
    if ~isscalar(d) || d <= 0
        error('d must be a positive scalar');
    end
    
    % 路径参数
    if ~isscalar(P) || P < 1 || mod(P,1) ~= 0
        error('P must be a positive integer');
    end
    
    % 用户参数
    if ~isscalar(K) || K < 1 || mod(K,1) ~= 0
        error('K must be a positive integer');
    end
    if length(angles) ~= K
        error('Angles structure must have K elements');
    end
    
    % 角度验证
    for k = 1:K
        if length(angles(k).gamma_x) ~= P || length(angles(k).gamma_y) ~= P
            error('Each user must have P gamma values');
        end
    end
end